####################################
# Main effect RDD
####################################

rm(list=ls())

library(gridExtra)
library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# reada data
d = read.dta13("~/Dropbox/Anti-american attitudes/01_data/clean_data/trump_election_data_clean_final.dta")   
names(d)

# countries
d_sal = d[d$pais==3,]
d_hon = d[d$pais==4,]
d_par = d[d$pais==12,]
d_ven = d[d$pais==16,]
d_dr = d[d$pais==21,]

############################
# US trust binary Salvador
############################

# bandwidths
h_us_trust_bin = rdbwselect(d_sal$us_trust_bin,d_sal$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d_sal$score, 0, bw = h_us_trust_bin, kernel = "triangular")
coef(summary(lm(d_sal$us_trust_bin ~ d_sal$treatment + d_sal$score + d_sal$treatment*d_sal$score, weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d_sal$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d_sal$us_trust_bin ~ d_sal$treatment + d_sal$score + d_sal$treatment*d_sal$score, weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d_sal$us_trust_bin ~ d_sal$treatment + d_sal$score + d_sal$treatment*d_sal$score, weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d_sal$us_trust_bin ~ d_sal$treatment + d_sal$score + d_sal$treatment*d_sal$score, weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d_sal$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="black"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)
results
mean(results$pe)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6)) + coord_cartesian(ylim=c(-0.6,0.6)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "black" = "black")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

############################
# US trust binary Honduras
############################

# bandwidths
h_us_trust_bin = rdbwselect(d_hon$us_trust_bin,d_hon$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d_hon$score, 0, bw = 7, kernel = "triangular")
coef(summary(lm(d_hon$us_trust_bin ~ d_hon$treatment + d_hon$score + d_hon$treatment*d_hon$score, weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d_hon$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d_hon$us_trust_bin ~ d_hon$treatment + d_hon$score + d_hon$treatment*d_hon$score, weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d_hon$us_trust_bin ~ d_hon$treatment + d_hon$score + d_hon$treatment*d_hon$score, weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d_hon$us_trust_bin ~ d_hon$treatment + d_hon$score + d_hon$treatment*d_hon$score, weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d_hon$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="black"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)
results
mean(results$pe)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6)) + coord_cartesian(ylim=c(-0.6,0.6)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "black" = "black")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

############################
# US trust binary Paraguay
############################

# bandwidths
h_us_trust_bin = rdbwselect(d_par$us_trust_bin,d_par$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d_par$score, 0, bw = h_us_trust_bin, kernel = "triangular")
coef(summary(lm(d_par$us_trust_bin ~ d_par$treatment + d_par$score + d_par$treatment*d_par$score, weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d_par$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d_par$us_trust_bin ~ d_par$treatment + d_par$score + d_par$treatment*d_par$score, weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d_par$us_trust_bin ~ d_par$treatment + d_par$score + d_par$treatment*d_par$score, weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d_par$us_trust_bin ~ d_par$treatment + d_par$score + d_par$treatment*d_par$score, weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d_par$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="black"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)
results
mean(results$pe)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6)) + coord_cartesian(ylim=c(-0.6,0.6)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "black" = "black")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

############################
# US trust binary Venezuela
############################

# bandwidths
h_us_trust_bin = rdbwselect(d_ven$us_trust_bin,d_ven$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d_ven$score, 0, bw = h_us_trust_bin, kernel = "triangular")
coef(summary(lm(d_ven$us_trust_bin ~ d_ven$treatment + d_ven$score + d_ven$treatment*d_ven$score, weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d_ven$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d_ven$us_trust_bin ~ d_ven$treatment + d_ven$score + d_ven$treatment*d_ven$score, weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d_ven$us_trust_bin ~ d_ven$treatment + d_ven$score + d_ven$treatment*d_ven$score, weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d_ven$us_trust_bin ~ d_ven$treatment + d_ven$score + d_ven$treatment*d_ven$score, weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d_ven$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="black"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)
results
mean(results$pe)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6)) + coord_cartesian(ylim=c(-0.6,0.6)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "black" = "black")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

#######################################
# US trust binary Dominican Republic
#######################################

# bandwidths
h_us_trust_bin = rdbwselect(d_dr$us_trust_bin,d_dr$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d_dr$score, 0, bw = h_us_trust_bin, kernel = "triangular")
coef(summary(lm(d_dr$us_trust_bin ~ d_dr$treatment + d_dr$score + d_dr$treatment*d_dr$score, weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d_dr$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d_dr$us_trust_bin ~ d_dr$treatment + d_dr$score + d_dr$treatment*d_dr$score, weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d_dr$us_trust_bin ~ d_dr$treatment + d_dr$score + d_dr$treatment*d_dr$score, weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d_dr$us_trust_bin ~ d_dr$treatment + d_dr$score + d_dr$treatment*d_dr$score, weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d_dr$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="black"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)
results
mean(results$pe)

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6)) + coord_cartesian(ylim=c(-0.6,0.6)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "black" = "black")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

